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Abstract 

Continuum  phase  field  theory  is  applied  to  study  elastic  twinning  in  calcite 
and  sapphire  single  crystals  subjected  to  indentation  loading  by  wedge-shaped 
indenters.  An  order  parameter  is  associated  with  the  magnitude  of  stress-free 
twinning  shear.  Geometrically  linear  and  nonlinear  theories  are  implemented 
and  compared,  the  latter  incorporating  neo-Hookean  elasticity.  Equilibrium 
configurations  of  deformed  and  twinned  crystals  are  attained  numerically 
via  direct  energy  minimization.  Results  are  in  qualitative  agreement  with 
experimental  observations:  a  long  thin  twin  forms  asymmetrically  under  one 
side  of  the  indenter,  the  tip  of  the  twin  is  sharp  and  the  length  of  the  twin  increases 
with  increasing  load.  Qualitatively  similar  results  are  obtained  using  isotropic 
and  anisotropic  elastic  constants,  though  the  difference  between  isotropic  and 
anisotropic  results  is  greater  in  sapphire  than  in  calcite.  Similar  results  are 
also  obtained  for  nanometer- scale  specimens  and  millimeter-scale  specimens. 
Indentation  forces  are  greater  in  the  nonlinear  model  than  the  linear  model 
because  of  the  increasing  tangent  bulk  modulus  with  increasing  pressure  in  the 
former.  Normalized  relationships  between  twin  length  and  indentation  force 
are  similar  for  linear  and  nonlinear  theories  at  both  nanometer  and  millimeter 
scales.  Twin  morphologies  are  similar  for  linear  and  nonlinear  theories  for 
indentation  with  a  90°  wedge.  However,  in  the  nonlinear  model,  indentation 
with  a  120°  wedge  produces  a  lamellar  twin  structure  between  the  indenter  and 
the  long  sharp  primary  twin.  This  complex  microstructure  is  not  predicted  by 
the  linear  theory. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Twinning  is  a  fundamental  deformation  mechanism  that  occurs  in  many  crystalline  solids.  In 
this  work,  phase  field  theory  is  used  to  model  deformation  twinning,  also  known  as  mechanical 
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twinning,  which  is  defined  as  twinning  induced  by  mechanical  stresses  [1-3].  Deformation 
twinning  is  henceforth  simply  referred  to  as  ‘twinning’.  The  transformation  strain  associated 
with  twinning  is  a  simple  shear.  Across  the  habit  plane,  the  orientation  of  the  Bravais  lattice 
differs  by  a  reflection  or  rotation  depending  on  the  kind  of  twin  under  consideration.  The 
sheared  and  reoriented  crystal  is  labeled  the  ‘twin’ ,  while  the  region  of  crystal  that  maintains  its 
original  lattice  orientation  is  labeled  the  ‘parent’  or  the  ‘matrix’ .  At  the  atomic  scale,  twinning 
takes  place  by  coordinated  movement  of  partial  dislocations  (i.e.  twinning  dislocations)  and/or 
atomic  shuffles  [2] .  A  finite  surface  energy  can  be  associated  with  the  twin  boundary,  often 
estimated  on  the  order  of  the  appropriate  stacking  fault  energy  [4] .  The  interfacial  energy  of  a 
growing  or  shrinking  twin  can  also  include  elastic  and  core  energies  of  twinning  dislocations 
comprising  such  interfaces,  often  labeled  ‘incoherent  interfaces’  [1].  Although  exceptions 
exist,  deformation  twinning  is  often  a  preferred  deformation  mechanism  in  crystals  with  low 
stacking  fault  energy  and/or  low  crystallographic  symmetry  [2] . 

Twinning  in  calcite  (CaCOs)  was  extensively  characterized  via  a  number  of  experiments  in 
the  mid-20th  century  [1, 5-11].  At  room  temperature,  dislocation  glide  does  not  occur  readily 
in  calcite,  but  twinning  can  occur  profusely  [  10] .  Thus,  calcite  is  an  ideal  material  for  validation 
of  models  of  deformation  twinning,  since  twin  morphologies  can  be  observed  using  optical 
measurements  (a  benefit  of  calcite ’s  transparency),  and  since  dislocation-mediated  plasticity 
need  not  be  addressed. 

A  number  of  experimental  studies  of  twinning  in  calcite  considered  indentation  loading, 
either  with  knife-edge  indenters  [5, 7, 9]  or  with  spherical  indenters  [8].  Calcite  crystals  were 
oriented  such  that  the  direction  of  twinning  shear  was  parallel  to  the  loading  direction  and 
the  habit  plane  normal  was  perpendicular  to  the  loading  direction.  Twinning  under  such 
conditions  was  often  reported  to  be  ‘elastic’,  that  is,  twins  under  the  indenter  disappeared 
fully  or  partially  upon  load  removal  [1, 5,  8, 11].  Twins  originated  at  the  contact  surface  and 
were  reported  to  be  thin  with  a  sharp  tip  or  apex  as  shown  in  figure  1,  with  the  length  of 
the  twin  increasing  with  increasing  load.  Twins  could  be  maintained  in  a  stable  position  (i.e. 
held  at  a  constant  length)  when  the  load  was  held  fixed,  with  measured  sizes  reported  on  the 
order  of  micrometers  to  hundreds  of  mm  depending  on  the  magnitude  of  applied  load  and 
type  of  indenter  geometry  [1, 5,  8, 11].  When  the  calcite  crystals  exhibit  no  surface  defects, 
twins  appear  only  after  a  threshold  load  is  applied,  whereas  twins  can  appear  immediately  (i.e. 
upon  application  of  a  minimal  load)  at  surface  imperfections  when  the  crystals  exhibit  such 
defects  [8].  The  title  of  this  work  is  motivated  by  the  transparency  of  calcite  (and  sapphire,  as 
discussed  later)  which  enables  direct  optical  observation  of  elastic  twinning  during  indentation 
loading.  Appearance  and  disappearance  of  elastic  twins  cannot  be  visually  observed  within 
opaque  materials  such  as  metals. 

Previous  analytical  models  have  offered  a  partial  description  of  physics  of  elastic  twinning 
under  concentrated  loading  but  cannot  be  categorized  as  fully  predictive  in  terms  of  twin 
morphology  or  elastic  fields.  Lifshitz  [12] — see  summary  of  model  in  English  language 
in  [1] — incorporated  a  nonlinear  elastic  law  (with  non-convex  energy  density)  and  showed 
that  a  two-dimensional  twin  of  finite  thickness  should  terminate  within  the  material  at  a  corner 
point  with  zero  aperture  angle,  i.e.  a  sharp  cusp.  This  theory  [1, 12]  does  not  incorporate 
dislocation  mechanics  or  any  quantitative  value  of  interfacial  energy  and  does  not  enable 
prediction  of  the  size  of  the  twin,  but  it  does  give  qualitative  information  regarding  the  twin’s 
equilibrium  shape. 

A  different  theory  developed  by  Kosevich  and  colleagues  [1, 9, 11]  assumes  a  priori  that 
the  twin  is  long  and  thin  and  can  be  modeled  as  a  one-dimensional  continuum  elasticity 
problem.  This  theory  considers,  for  equilibrium,  a  balance  of  forces  of  linear  elastic  origin, 
Peierls  forces  [13]  associated  with  lattice  friction,  and  surface  tension  forces  that  resist  the 
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indentation  force  P 


Figure  1.  Wedge-shaped  twin  in  calcite,  after  Kosevich  and  Boiko  [1].  Twin  appears  only  under 
left  side  of  indenter  because  shear  deformation  under  right  side  is  directed  in  anti-twinning  sense. 


formation  of  stacking  faults  in  the  region  of  extension  of  the  twin  (i.e.  at  the  tip).  Linear  elastic 
forces  are  assumed  to  arise  from  a  continuous  distribution  of  twinning  (partial)  dislocations 
comprising  the  interface.  The  thickness  of  a  one-dimensional  twin  at  a  point  is  related  to  the 
dislocation  density  distribution.  Phenomenological  parameters  in  this  theory  can  be  prescribed 
so  that  the  analytical  solution  matches  experimental  measurements  of  twin  length  and  applied 
indentation  force  [1, 9, 1 1],  but  the  theory  cannot  be  regarded  as  fully  predictive  because  such 
parameters  are  not  obtained  from  independent  measurements. 

Crystal  plasticity-type  models  have  also  been  developed  to  describe  twinning.  In  such 
models,  the  volume  or  mass  fraction  of  a  given  material  element  occupied  by  twins  belonging 
to  one  or  more  twin  systems  is  treated  as  an  internal  state  variable  that  evolves  with 
thermomechanical  loading  via  a  kinetic  law.  Transformation  kinetics  are  typically  controlled 
by  driving  stress  and/or  temperature.  These  models  are  useful  for  describing  texture  evolution 
[14, 15]  and  macroscopic  stress-strain  behavior  [3, 16-18].  A  crystal  mechanics  approach  has 
also  been  used  to  predict  subgrain  twin  morphology  in  3D  poly  crystal  simulations  [19]. 

In  this  paper,  a  phase  field  theory  for  twinning  developed  recently  by  the  authors  [20]  is 
used.  In  any  given  simulation,  twinning  is  restricted  to  occur  on  only  one  twin  system.  An  order 
parameter  gauges  the  magnitude  twinning  shear  on  this  system  at  any  location  in  the  problem 
domain.  The  general  form  of  the  model  in  [20]  accounts  for  large  deformations,  nonlinear 
and  anisotropic  mechanical  elastic  properties,  and  anisotropic  surface  energy  associated  with 
twin  boundaries.  Equilibrium  configurations  of  deformed  and  twinned  crystals  are  attained 
via  direct  energy  minimization.  The  theory  is  framed  in  the  null  temperature  (i.e.  OK)  limit, 
with  no  account  of  dissipation  or  irreversibility.  Model  assumptions  of  static  equilibrium 
and  null  dissipation  are  deemed  particularly  appropriate  for  a  study  of  ‘elastic’  twinning  in 
indented  calcite.  Apart  from  the  geometry  of  the  twin  system  (shearing  plane,  direction  and 
magnitude)  known  a  priori  from  the  crystal  structure,  the  only  material  properties  entering  the 
present  phase  field  model  are  the  elastic  constants  which  are  known  from  experiments  [21], 
and  the  twin  boundary  energy  and  characteristic  thickness,  which  are  known  from  atomic 
calculations  [22].  Hence,  the  model  is  deemed  fully  predictive  because  morphologies  and 
elastic  fields  associated  with  twinning  do  not  result  from  any  a  priori  assumptions  on  twin 
shape  or  kinetics  of  inelastic  material  behavior. 

In  phase  field  modeling,  interfaces  between  pure  phases — in  the  present  case  between 
twinned  and  untwinned  regions  of  crystal — possess  a  finite  thickness.  Boundary  value 
problems  incorporating  diffuse  interface  theory  are  amenable  to  treatment  with  modem 
numerical  methods  [23,24].  In  contrast,  sharp  interface  models  that  address  equilibrium  or 
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stability  among  phases  of  nonlinear  elastic  materials  (including  deformation  twins)  have  been 
examined  extensively  using  analytical  methods  [25-29],  and  less  so  with  numerical  techniques. 

In  addition  to  calcite,  sapphire  is  also  considered.  Sapphire  is  the  single  crystal  form  of 
a -alumina  (AI2O3),  and  is  alternatively  referred  to  in  various  literature  as  corundum.  High- 
purity  sapphire,  like  calcite,  is  transparent.  Results  for  sapphire  provide  a  useful  basis  for 
comparison  with  those  for  calcite.  Sapphire  is  notably  stiffer  (much  higher  hardness  and 
elastic  constants)  than  calcite.  Calcite  and  sapphire  belong  to  the  same  crystallographic  space 
group;  both  are  trigonal.  Twinning  in  sapphire  can  occur  on  either  of  two  twin  systems  (on  basal 
or  rhombohedral  planes)  [16, 18,  30, 31];  both  twin  systems  are  considered  individually  in  this 
work.  Dislocation  plasticity  (i.e.  plastic  slip)  can  accompany  twinning  in  sapphire  subjected  to 
indentation  loading  [32-34],  since  critical  resolved  shear  stresses  required  to  induce  twinning 
and  slip  are  of  the  same  order  of  magnitude  [16, 18].  Because  the  model  and  results  presented 
here  do  not  consider  plastic  slip,  they  represent  idealized  behavior.  The  present  model  would 
require  substantial  changes  to  address  kinematics  and  kinetics  of  dislocation  distributions  and 
plastic  slip  [3, 16-18,  35],  with  corresponding  dissipation. 

Reported  first  in  this  paper  are  predictions  for  linear  elastic  behavior.  In  the  present 
phase  field  model  [20],  when  elastic  strains  are  small,  the  difference  in  driving  shear  stress  for 
twinning  is  a  factor  on  the  order  of  twinning  shear.  Previous  theoretical  studies  of  homogeneous 
twin  nucleation  [2,36,37]  have  relied  on  linear  elasticity,  presumably  because  analytical 
solutions  are  not  feasible  for  nonlinear  models.  The  geometrically  linear  model  considered 
in  this  work  enables  direct  comparison  with  analytical  solutions  of  wedge  indentation  in 
elastic  and  elastic-plastic  materials  [38-40] .  Geometrically  linear  theory  has  also  been  applied 
elsewhere  towards  phase  field  modeling  of  twinning  [41]. 

Reported  second  in  this  paper  are  model  predictions  for  a  finite  strain  model  incorporating 
nonlinear  elastic  behavior.  Results  from  nonlinear  and  linear  theories  are  also  compared 
directly.  Linear  theory  can  provide  qualitative  insight  and  in  some  cases  quantitatively 
reasonable  descriptions  of  phase  transformations  and  twinning  behavior  [29].  However,  in 
other  cases,  linear  and  nonlinear  models  can  yield  drastically  different  results  with  regards  to 
both  analytical/theoretical  predictions  [29]  and  numerical  predictions  of  phase  field  models 
[42].  Results  reported  in  this  paper  demonstrate  qualitative  agreement  between  linear  and 
nonlinear  models  for  some  phenomena,  but  significant  differences  for  others.  The  theory 
and  results  presented  here  incorporate  a  neo-Hookean  model  [43, 44]  for  nonlinear  elastic 
behavior,  thus  presuming  an  isotropic  elastic  response.  Such  an  assumption  is  thought  justified 
since  differences  in  predictions  among  isotropic  and  anisotropic  linear  elastic  models  are  not 
drastic,  with  qualitatively  similar  twin  sizes  and  shapes  obtained  for  isotropic  and  anisotropic 
elasticity. 

The  rest  of  this  paper  is  organized  as  follows.  The  phase  field  theory  is  summarized  in 
section  2,  applicable  to  any  generic  crystalline  material  undergoing  the  same  deformation 
mechanisms.  Properties  of  calcite  and  sapphire  are  addressed  in  section  3.  Numerical 
application  of  the  linear  theory  to  indentation  is  discussed  in  section  4;  application  of  the 
nonlinear  theory  in  section  5.  Conclusions  follow  in  section  6. 

Notation  of  continuum  mechanics  is  used.  Familiarity  of  the  reader  with  elasticity  theory 
and  phase  field  modeling  is  assumed;  overviews  of  the  latter  include  [23,24].  Vectors 
and  higher-order  tensors  are  written  in  bold  italic  font;  scalars  and  components  of  vectors 
and  tensors  are  written  in  italic  font.  When  indicial  notation  is  used,  summation  proceeds 
over  repeated  indices.  Vectors  and  tensors  are  referred  to  a  fixed  Cartesian  frame  of 
reference,  with  indices  in  the  subscript  position.  The  scalar  product  of  vectors  a  and  b 
is  a  b  =  apJoA  —  a\b\  +  a^b^  +  03^3  in  three-dimensional  space.  The  outer  product  is 
(a  0  b)AB  =  cla^b-  Juxtaposition  implies  summation  over  one  set  of  adjacent  indices, 
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e.g.,  ( AB)ab  =  AacBCb •  The  colon  denotes  summation  over  two  sets  of  indices,  e.g., 
A  :  B  =  AabBab •  A  T  superscript  denotes  the  transpose,  e.g.,  Atab  =  ABa- 


2.  Theory 

A  continuum  phase  field  theory  for  twinning  [20]  is  summarized.  The  theory  for  geometrically 
linear  kinematics  is  presented  first,  followed  by  the  finite  strain  theory. 


2.7.  Linear  theory 


Let  X  g  £2  be  a  material  point.  The  order  parameter  function  is  where  t  denotes  time. 

The  order  parameter  distinguishes  between  two  phases:  (i)  the  original  crystal  (the  parent)  and 
(ii)  the  twin.  Interfaces  between  phases  are  twin  boundaries.  Order  parameter  f)  generally 
exhibits  the  following  values: 

rj(X,  •)  =  0VX  e  parent, 

=  1VX  g  twin, 

e  (0,  1)VX  g  twin  boundaries.  (1) 


In  linear  elasticity,  kinematic  field  variables  are  displacement  u  and  its  gradient 

f3  =  Vm,  Pab  =  Vbua.  (2) 

The  distortion  is  split  additively  as 

/3  =  /3e  +  /3",  (3) 

where  /3E  is  the  elastic  distortion  and  (3,]  is  the  stress-free  distortion  associated  with  the  twinning 
shear  or  twinning  transformation: 

/3''  =  [<?(/?)]  y0  s  ®  m.  (4) 

The  unit  normal  to  the  surface  of  composition  (i.e.  the  habit  plane)  is  m.  The  magnitude  of 
the  twinning  shear  and  the  shear  direction  are  yo  and  s,  respectively.  The  following  conditions 
hold:  m  -  m  —  s  s  —  1  and  s  •  m  =  0. 

Interpolation  function  ip  relates  twinning  distortion  and  the  order  parameter  [42]: 

<p{ri)  =  342  -  2^3.  (5) 

Note  that  ip  is  a  monotonically  increasing  function  obeying  ip{ 0)  =  0,  ip{\)  =  1 ,  with  vanishing 
derivatives  with  respect  to  r\  at  its  endpoints:  <^(0)  =  ip' (l)  =  0.  This  function  also  obeys  the 
anti-symmetry  conditions  ip (l  —  r])  =  1  —  (p(rj).  The  symmetric  elastic  strain  tensor  is 

eE(Vu,  n)  =  i(/3E  +  /3ET) 

=  \  { Vu  +  (Vw)T  —  yo[(p(rj)][s  <g>  m  +  m  <g>  s]}  .  (6) 

The  total  free  energy  functional  for  a  body  of  volume  Q  is  written  as 

^(u,r])=  I  W(Vu,  r])  d£2  +  /  / (rj ,  V r])  d£2 ,  (7) 

J  £2  J  £2 

where  W  is  the  elastic  strain  energy  density  and  /  accounts  for  interfacial  energy.  Strain 
energy  density  and  second-order  elastic  moduli  are 

w  =  W[eE(Vu,  n),  /?]  =  yEABCABCD(r])eED,  (8) 


CAbcd(i1 )  = 


d2W 


eE=0 


(9) 
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For  a  centro symmetric  structure,  a  re-orientation  matrix  Q  associated  with  twinning  is  the 
reflection  [2, 3] 


Q  =  1  —  2 m  (g)  ra,  (10) 

with  1  the  unit  tensor.  Elastic  coefficients  of  the  fully  twinned  crystal  are  related  to  those  of 
the  parent  by 

Cabcd(  1)  =  QaeQbfQcgQdhCefgh(Q )•  (11) 

Elastic  coefficients  in  interfacial  regions  are  interpolated  as 

Cabcd (xi)  =  Cabcd( 0)  +  [Cabcd(  1)  —  Cabcd( 0)]  (p(rf).  (12) 

When  isotropic  elasticity  is  used,  elastic  coefficients  are  independent  of  orientation: 

Cabcd  =  ^ab^cd  +  ^(8ac$bd  +  8ad$bc),  (13) 

with  X  and  [i  the  Lame  constant  and  shear  modulus;  hence  in  the  isotropic  case  W  does  not 
explicitly  depend  on  r]. 

The  local  interfacial  energy  per  unit  volume  is 

f(n,  Vri)  =  f0(rj)  +  n(rj)  :  (Wrj  0  V*j),  (14) 

with  hi  a  symmetric  second-order  tensor.  For  anisotropic  surface  energy,  kab(J])  =  «ab( 0)  + 
$(p(rj)\KAB(\)  ~  kab( 0)],  where  similarly  to  (12),  tcAB(  1)  =  QacQbdKcd( 0),  and  where 
f  G  [0,  1]  is  a  scalar  constant.  Case  f  =  0  considered  in  previous  work  [20]  implies  that  the 
gradient  contribution  to  interfacial  energy  is  invariant  with  respect  to  rj  and  hence  the  twinning 
shear;  condition  f  =  1  implies  that  the  gradient  contribution  depends  on  the  transformation 
of  the  lattice  (e.g.  reflection)  associated  with  twinning.  Other  models  have  included  highly 
anisotropic  surface  energy  [41].  When  interfacial  energy  is  isotropic,  k  =  tel  and 

/0j>  Vtj)  =  /o(?7)  +  k|Vt7|2.  (15) 

Prescribed  for  /0  is  a  standard  ‘double- well’  potential: 

/o(»?)  =  V(  1  -  >7)2,  (16) 

with  A  >  0.  In  the  isotropic  approximation  A  and  /c  are  related  to  equilibrium  energy  per  unit 
area  T  and  thickness  l  of  an  unstressed  interface  via  [20] 

k  =  3T//4,  A  =  12E/Z.  (17) 

In  this  work,  both  /0  and  V77  contribute  to  what  is  called  interfacial  energy  (14).  This  label 
follows  the  scheme  denoted  in  (1):  by  the  convention  used  here,  values  of  77  ^  0,  1  designate 
‘interface’,  so  /  =  0  wherever  r\  =  0,1.  Other  presentations  may  use  different  notation, 
e.g.  fo  has  been  labeled  ‘deformation  energy’  in  another  phase  field  model  for  twinning  [41], 
since  a  uniform  value  of  the  order  parameter  differing  from  zero  or  unity  can  result  in  a  finite 
energetic  contribution.  It  is  also  noted  that  (15)  with  (16)  constitute  a  simple  assumption 
regarding  interfacial  energy;  more  elaborate  models  accounting  for  additional  physics  (e.g. 
more  than  two  variants,  anisotropic  dislocation  core  and  stacking  fault  contributions,  or  grain 
boundary  misorientation  and  cusps  associated  with  low-energy  boundaries)  exist  [41, 45-48]. 

The  following  variational  equation  is  posited  that  will  yield  equilibrium  equations  and 
boundary  conditions: 

8^  =  <b  t-8udS+(b  h8r]dS,  (18) 

JdQ  JdQ 
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where  t  is  a  mechanical  traction  vector  per  unit  area,  d S  is  a  surface  element  of  3  £2  and  h  is  a 
scalar  conjugate  force  to  variations  of  the  order  parameter.  Euler-Lagrange  equations  obtained 
from  (18)  are  [20] 


V  • 


3  W 
3  Vix 


=  V  •  <7  =  0, 


/<$ 


dW 

2V  •  (kWt})  + -  =  0, 

dr]  Vu 


(19) 


where  a  is  the  symmetric  stress  tensor.  Corresponding  boundary  conditions  are 


t  =  an,  h  =  2k,  :  (V77  (g)  n), 


(20) 


where  n  is  the  unit  outward  normal  to  3  £2.  The  stress  also  obeys  the  usual  linear  elasticity 
relationship 

&ab  =  Cabcd  ^cd-  (21) 

The  second  of  equilibrium  conditions  (19)  can  be  rewritten  as 
fo  =  2Va(kab^bti)  +  <p'{rYo  +  ^aV^bV  [kAb( 0)  —  kab(  1)1 

+  (\/2)eEABACD  [Cabcd( 0)  -  CABCD(m,  (22) 

where  the  driving  shear  stress  for  twinning  is 


r  =  a  :  (s  0  m). 


(23) 


For  isotropic  surface  and  elastic  energies,  with  <p  from  (5)  and  /o  from  (16),  equilibrium 
condition  (22)  becomes 

3r/o^(l  —  rj)  =  At]  (l  —  3  rj  +  2  r]2)  —  kW2t].  (24) 

Both  sides  of  (22)  and  (24)  vanish  in  regions  of  uniform  phase  where  rj  =  0  or  rj  =  1. 


2.2.  Nonlinear  theory 

Equations  (1),  (5)  and  (14)— (17)  also  apply  for  the  nonlinear  theory.  Spatial  x  and  referential 
X  coordinates  of  a  material  particle  are  related  by  the  motion 

x  =  X(X,t).  (25) 

The  deformation  gradient  is 

F  =  V*.  FaA  =  VAXa,  (26) 

with  VA  =  d/dXA.  The  deformation  gradient  is  decomposed  multiplicatively  as 

F  =  FeF\  (27) 

where  FE  =  FF'r'  is  the  elastic  deformation  and  F^lrjiX,  t)  \  represents  the  stress-free 
twinning  shear.  For  twinning  on  a  single  system,  the  following  form  applies: 

Fv  =  1  +  [<p(rj) ]  y0  s®m.  (28) 

The  unit  (identity)  tensor  is  1.  The  symmetric  elastic  deformation  tensor  is 

Ce  =  FeTFe,  CEp  =  FEaFEp.  (29) 

Twinning  is  isochoric: 

J71  =  detF^  =  1  +  [<p(ri)]  y0  s  -  m  =  1.  (30) 

It  follows  that 

J  =  det  F  =  det  FE  det  Fn  =  JET  =  JE  =  Vdet  CE.  (31) 
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The  total  energy  functional  for  a  body  of  initial  (i.e.  undeformed)  volume  Q  is 

'I'(X,77)=  /  W(F,  rj)dQ+  f  f(rj,Vrj)  dQ,  (32) 

Jn  Jo. 

where  W  is  the  elastic  strain  energy  density  and  /  accounts  for  interfacial  energy.  For  a 
neo-Hookean  material  [44]  ,  strain  energy  density 

W  =  W[Ce(F,  rj)]  =  ^(tr  CE  -  3)  +  In  J  Q  In/  -  ]ij  .  (33) 

For  a  (deviatoric)  simple  shear  deformation  gradient,  the  response  of  a  neo-Hookean  material 
is  the  same  as  that  for  a  linear  elastic  material  with  the  same  shear  modulus.  On  the  other  hand, 
when  large  volume  changes  are  involved,  predictions  of  neo-Hookean  and  linear  elastic  models 
can  differ  significantly.  Specifically,  as  will  be  demonstrated  later  in  section  3,  a  neo-Hookean 
model  can  describe  the  usual  increase  in  tangent  bulk  modulus  with  increasing  compressive 
pressure. 

The  following  variational  equation  is  posited  that  will  yield  equilibrium  equations  and 
boundary  conditions: 


8V  =  (f  t-8XdS+(fi  h  8rj  dS,  (34) 

where  t  is  a  mechanical  traction  vector  per  unit  reference  area,  d S  is  a  surface  element  of 
3  £2,  and  h  is  a  scalar  conjugate  force  to  variations  of  the  order  parameter.  Application  of  the 
divergence  theorem,  integration  by  parts,  and  localization  of  (34)  lead  to  the  Euler-Lagrange 
equations  [20]: 


dw 

Jf 


=  v  •  p  =  o, 


/o-2kVV 


dw 

3  rj 


=  0, 


(35) 


with  P  the  first  Piola-Kirchhoff  stress.  Corresponding  boundary  conditions  are 


t  =  Pn ,  h  =  2 kVt]  •  n, 


(36) 


where  n  is  the  unit  outward  normal  to  3  Q. 
The  first  Piola-Kirchhoff  stress  satisfies 


3  W 

Jf 


3  W  3  Ch  3  Fh 


3  CE  '  3  FE  '  3  F 


-T 


The  following  kinematic  identities  prove  useful  for  neo-Hookean  elasticity: 


d  _  l^E-1  ax  _  _  E_i  E_i  ^  E-i  E-i 

3 CE  ~  2  9  3 “  ab  Px  xb  afi' 

The  symmetric  elastic  second  Piola-Kirchhoff  stress  is  obtained  from  (33)  and  (37)  as 


E  =  2—^  =jul  +  (UnJ-/x)CE  . 

The  second  of  equilibrium  conditions  (19)  can  be  rewritten  as 
fo  =  2/cV2r]  +  ryo<p', 


where  the  driving  shear  stress  for  twinning  is 


1  3  W 


x  = 


Yo<P'  dr] 
=  E  :  [CE(s 


1  dW  dCE(F,  i]) 


f  y0<p'  dCE 


dr] 


(37) 


(38) 


(39) 

(40) 


Using  (5),  (16)  and  (41),  equilibrium  condition  (40)  can  then  be  written  as  in  (24). 


(41) 
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It  can  be  shown  [20],  given  certain  boundary  conditions,  that  a  (meta)stable  configuration 
of  a  body  undergoing  twinning  corresponds  to  a  (local)  minimizer  of  energy  functional 
of  (32).  Thus  the  mathematical  problem  of  interest  whose  solution  is  sought  numerically,  as 
described  later  in  sections  4  and  5,  can  be  stated  simply  as 

min^(x,  rj).  (42) 


3.  Materials 

Relevant  behaviors  of  calcite  and  sapphire  are  discussed;  requisite  material  properties 
associated  with  elasticity  and  twinning  in  each  material  are  tabulated. 


3.1.  Calcite 


Calcite  is  a  stable  low-temperature  polymorph  of  calcium  carbonate  (CaCOs),  and  is  an  ionic 
crystal  with  trigonal  (i.e.  rhombohedral)  symmetry  [49,50].  Transparent  single  crystals  are 
also  called  iceland  spar.  Information  regarding  the  primitive  unit  cell  and  requisite  properties  of 
calcite  are  listed  in  table  1  with  supporting  references.  Calcite  belongs  to  space  group  R3 c  and 
centrosymmetric  point  group  3m.  The  primitive  rhombohedral  unit  cell  is  described  by  lattice 
parameter  a  and  bond  angle  a.  The  primitive  cell  contains  ten  atoms  (two  CaCC>3  formula 
units);  several  additional  parameters  are  needed  to  fully  describe  the  polyatomic  structure 
and  can  be  found  elsewhere  [50].  Hexagonal  crystallographic  notation  is  sometimes  used  to 
describe  calcite  [14,  36, 50],  which  exhibits  four  different  kinds  of  twin  systems  [22].  For  the 
e+  system  considered  here,  the  shear  is  large:  yo  =  0.694.  Miller  indices  shown  in  figure  1 
for  shearing  direction  and  habit  plane,  (00  1)  { 1  10},  correspond  to  the  cleavage  rhombohedral 
pseudocell  [51]  rather  than  the  primitive  unit  cell  of  the  Bravais  lattice. 

Twin  boundary  surface  energy  is  obtained  from  lattice  statics  (0  K)  calculations  conducted 
elsewhere  [22]  using  an  empirical  potential  [52].  In  most  phase  field  simulations  discussed 
later,  surface  energy  is  assumed  isotropic,  following  usual  theoretical  studies  [36,  37],  but  in 
one  case  anisotropy  of  the  surface  energy  is  considered.  In  a  material  coordinate  system  with 
axes  aligned  parallel  to  twinning  direction  and  habit  plane  normal  (Xi||m  and  X2II s),  the 
gradient  coefficient  entering  (14)  is  written  [20, 41]  as 


Kn  0 
0  /C22 


(43) 


In  the  isotropic  case, /ci  1  =  K22  =  k.  The  anisotropic  case  considered  later  is  K22  =  4/cn  =  4/c, 
which  would  seem  to  favor  boundaries  extended  parallel  to  the  habit  plane.  Recall  from  (17) 
that  surface  energy  T  a  /c1/2.  The  rationale  for  anisotropic  surface  energy  is  that  twinning 
dislocations  at  a  moving  portion  of  the  boundary  (i.e.  an  incoherent  interface  [1, 20, 41])  would 
contribute  core  and  elastic  energies  to  the  total  surface  energy  of  the  interface.  In  contrast,  the 
fully  formed  (i.e.  coherent)  twin  boundary  surface  would  have  less  energy  than  such  a  moving 
portion  because  it  does  not  contain  energy  of  dislocations,  just  stacking  fault  energy  associated 
with  reflection  of  the  lattice  across  the  twin  boundary.  Anisotropic  values  considered  here  are 
anticipated  to  cause  the  twin  to  elongate  in  the  direction  of  s  and  shorten  in  the  direction  of 
m,  in  order  to  decrease  the  contribution  of  /C22(V2?7)2  to  the  energy  in  (14).  It  is  also  assumed 
that  hi  does  not  depend  on  r]  and  hence  f  =  0  in  (22). 

Equilibrium  thickness  l  over  which  atoms  deviate  from  their  ideal  positions  is  chosen  as 
1  nm,  corresponding  to  about  five  {110}  planes.  This  value  follows  from  atomic  simulations 
of  twin  boundaries  in  calcite  [22],  where  it  is  reported  that  the  relaxed  crystal  structure  differs 
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Table  1.  Properties  of  calcite  single  crystals. 


Parameter 

Value 

Definition 

Reference 

a 

6.201 

lattice  parameter  (A) 

[50] 

a 

48.3° 

rhombohedral  angle 

[50] 

Yo 

0.694 

shear  for  e+  twin 

[14, 36] 

C 11 

165.4(148.1) 

elastic  constants  0  K  (300  K)  (GPa) 

[21] 

C12 

65.0(55.8) 

C13 

61.6(54.6) 

C14 

— 22.8(— 20.6) 

C33 

89.5(85.6) 

C44 

36.6(32.7) 

X 

61.7(54.6) 

Lame  modulus  0  K  (300  K)  (GPa) 

40.2(36.7) 

shear  modulus  0  K  (300  K)  (GPa) 

K 

88.5(79.1) 

bulk  modulus  0  K  (300  K)  (GPa) 

(73.5) 

[54] 

K' 

4.0 

(dK/dp)p= 0 

[54] 

r 

183 

twin  boundary  energy  (mJ  m-2) 

[22] 

/ 

1.0 

twin  boundary  thickness  (nm) 

[22,42] 

from  that  of  the  bulk  mainly  over  five  atomic  layers.  The  same  characteristic  thickness  value 
(1  nm)  has  also  been  used  in  phase  field  models  of  other  crystalline  materials  [20, 42]. 

Voigt’s  notation  [3, 53]  is  used  for  elastic  constants  in  table  1.  In  simulations  discussed 
later,  both  isotropic  and  anisotropic  elasticity  models  are  considered.  In  the  former,  Voigt- 
averaged  elastic  constants  [3, 4, 53]  are  used: 

*  =  ^(C11+5C12  +  8C13  +  C33-4C44),  (44) 

M  =  ^j(7Cn  —  5Ci2  —  4Ci3  +  2C33  +  I2C44).  (45) 

Elastic  anisotropy  of  calcite  is  significant:  C33  ~  0.5Cn,  C\\  ~  — O.6C44  and  C44  ~  0.7 C^. 
Recall  also  that  for  the  trigonal  crystal  class  of  calcite,  C22  =  C n,  C24  =  —  Cm  and 
2C*66  =  C\\  —  C 12.  In  simulations  that  follow,  elastic  constants  measured  at  160  K  [21] 
are  extrapolated  to  0  K  for  use  in  idealized  simulations  of  twinning  at  the  nanoscale  (indenter 
size  on  the  order  of  nm),  in  strict  adherence  to  the  theoretical  model  of  section  2  and  [20] 
that  does  not  consider  dissipation  or  thermal  phenomena.  Zero  temperature  solutions  are  also 
more  amenable  to  comparison  with  predictions  of  lattice  statics  [22] .  On  the  other  hand,  room 
temperature  elastic  constants  are  used  in  later  simulations  of  indentation  at  the  laboratory  scale 
(indenter  size  on  the  order  of  mm). 

To  validate  the  use  of  a  neo-Hookean  model  for  calcite,  the  response  of  the  material 
under  purely  volumetric  deformation,  F  =  /1/31,  is  computed  and  compared  with  the 
linear  elastic  response  and  that  for  the  Birch  equation  of  state  (EOS)  [55],  where  the  latter 
incorporates  the  bulk  modulus  K  and  the  pressure  derivative  of  the  bulk  modulus  K'  measured 
experimentally  [54].  Specifically,  Cauchy  pressure  p  is  computed  for  each  elastic  constitutive 


model  via 

p  =  K(l-J) 

(linear  elasticity) , 

(46) 

p  =  -  XJ~l  InJ 

(neo-Hookean) , 

(47) 

p  =  ^/“7/3(  1  -  y2/3)  1  +  ^(K'  -  4 )(J~2/3  -  1) 

(Birch  EOS). 

(48) 

Results  of  the  calculations  for  calcite  are  compared  in  figure  2(a).  In  each  case,  the  pressure  is 
normalized  by  the  same  value  of  Voigt- averaged  bulk  modulus  from  table  1,  Ky  =  79.1  GPa. 
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a 


(a) 


J 


(b) 


J 


Figure  2.  Cauchy  pressure  (normalized  by  Voigt-averaged  bulk  modulus)  versus  volume  for  linear 
and  nonlinear  elasticity  models  of  (a)  calcite  and  ( b )  sapphire. 


Close  agreement  is  observed  between  predictions  from  the  neo-Hookean  model  to  be  used  in 
simulations  and  those  from  the  Birch  EOS  that  is  fit  to  experimental  data  [54].  Note  that  the 
response  for  each  of  the  two  nonlinear  models  is  stiffer  than  that  of  linear  elasticity  at  higher 
pressures. 

A  polymorphic  phase  transition  from  trigonal  to  a  monoclinic  structure  has  been  observed 
in  calcite  deformed  under  hydrostatic  compression  at  pressures  on  the  order  of  1 .5  GPa  [54, 56] . 
Evidence  of  such  a  phase  transformation  has  not  been  reported  for  indentation  experiments 
on  calcite  [1,5, 7-11].  However,  it  is  plausible  that  such  a  transformation  could  occur,  given 
the  high  pressures  that  exist  immediately  beneath  the  indenter  that  can  easily  exceed  the 
transformation  pressure.  For  a  wedge  indenter  with  a  sharp  tip,  the  analytical  elastic  solution 
(linear  elasticity,  in  the  absence  of  twinning)  suggests  that  the  pressure  immediately  beneath 
the  tip  of  the  wedge  tends  to  infinity  [38, 40].  The  phase  field  theory  and  simulations  described 
in  this  work  do  not  address  the  possibility  of  this  phase  transformation  in  calcite. 


3.2.  Sapphire 

Sapphire  is  a  stable  low-temperature  polymorph  of  alumina  (AI2O3),  and  like  calcite  exhibits 
trigonal  symmetry  [18,30].  Sapphire  is  alternatively  known  as  a-A^C^  or  corundum; 
nominally  colorless  when  pure,  sapphire  with  trace  amounts  of  secondary  elements  often 
exhibits  blue  or  red  tints,  with  red  corundum  commonly  called  ruby.  Information  regarding 
the  primitive  unit  cell  and  requisite  properties  of  sapphire  are  listed  in  table  2.  Like  calcite, 
sapphire  belongs  to  space  group  R3c  and  centrosymmetric  point  group  3m.  The  primitive  cell 
contains  ten  atoms  (two  AI2O3  formula  units);  several  additional  parameters  are  needed  to  fully 
describe  the  polyatomic  structure  [30] .  Hexagonal  crystallographic  notation  is  typically  used 
to  describe  sapphire  [16, 18,30],  which  exhibits  two  kinds  of  twin  systems  [16, 18,30,31]. 
The  basal  (B)  system  corresponds  to  (1  I  00)  {000  1}  in  the  structural  unit  cell  [30],  and  the 
shear  is  large:  yo  =  0.635.  The  rhombohedral  (R)  system  corresponds  to  (1 1  0  T) { 1  I  02}  in 
the  structural  cell,  and  the  shear  is  somewhat  smaller:  yo  =  0.202. 

Again,  both  isotropic  and  anisotropic  elasticity  models  are  considered.  In  the  former, 
Voigt- averaged  elastic  constants  [53]  are  used.  Elastic  anisotropy  of  sapphire  is  mild: 
C33  ~  C 11,  C 14  ~  O.I6C44  and  C44  ~  0.9C66-  In  the  simulations  that  follow,  temperature 
dependence  of  elastic  coefficients  where  available  [57]  is  used  to  provide  0  K  constants  for 
use  at  the  nanoscale;  room  temperature  constants  [58]  are  used  in  simulations  of  indentation 
at  the  laboratory  (mm)  scale.  Following  [3, 4, 16, 18],  surface  energies  of  twin  boundaries  are 
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Table  2.  Properties  of  sapphire  single  crystals. 


Parameter 

Value 

Definition 

Reference 

a 

5.12 

lattice  parameter  (A) 

[30] 

a 

55.3° 

rhombohedral  angle 

[30] 

Yo 

0.635 

shear  for  basal  (B)  twin 

[30] 

0.202 

shear  for  rhombohedral  (R)  twin 

[31] 

Cn 

500.0  (498.0) 

elastic  constants  OK  (300K)  (GPa) 

[57,58] 

C 12 

167.8(163.0) 

C 13 

120.5(117.0) 

C\4 

23.7  (23.0) 

C33 

502.0  (502.0) 

C44 

151.0(147.0) 

X 

146.7  (144.2) 

Lame  modulus  0  K  (300  K)  (GPa) 

u 

166.5(165.7) 

shear  modulus  0  K  (300  K)  (GPa) 

K 

257.7  (254.7) 

bulk  modulus  0  K  (300  K)  (GPa) 

(254.4) 

[61] 

(239.0) 

[62] 

K' 

4.3 

(dK/dp)p=0 

[61] 

0.9 

[62] 

r 

745 

B  twin  boundary  energy  (mJ  m-2) 

[60] 

125 

R  twin  boundary  energy  (mJ  m-2) 

[59] 

/ 

1.0 

twin  boundary  thickness  (nm) 

[63,42] 

estimated  as  half  of  corresponding  stacking  fault  energies,  the  latter  determined  experimentally 
[59]  or  through  atomic  modeling  [60]. 

Use  of  a  neo-Hookean  model  for  sapphire  is  validated  for  volumetric  compression 
deformation  in  figure  2(b),  where  pressures  are  computed  as  in  section  3.1.  Pressures  obtained 
from  the  neo-Hookean  model  are  intermediate  to  those  from  the  Birch  EOS  fitted  to  two 
different  sets  of  hydrostatic  compression  experiments  [61, 62].  The  Birch  EOS  with  K  and  K' 
from  [62]  gives  predictions  very  similar  to  a  linear  elastic  model  incorporating  Voigt- averaged 
bulk  modulus  Ky  =  254.7  GPa. 

Although  inconclusive  evidence  exists  from  shock  compression  experiments  [64],  sapphire 
may  undergo  a  solid-solid  phase  transformation  at  pressures  in  excess  of  79  GPa.  Evidence 
of  such  a  transformation  has  not  been  reported  in  indentation  experiments  involving  twinning 
[32-34] ,  though  it  is  conceivable  for  such  high  pressures  to  be  achieved  for  very  sharp  indenters. 
The  possibility  of  such  a  high-pressure  phase  transformation  in  sapphire  is  omitted  in  the  present 
phase  field  simulations. 


4.  Phase  field  simulations:  linear  theory 

Simulations  of  wedge  indentation  into  calcite  and  sapphire  using  the  linear  theory  of  section  2. 1 
are  described.  Numerical  results  are  analyzed. 

4.1.  Boundary  value  problem 

The  problem  of  study  is  illustrated  in  figure  3.  Simulations  are  two-dimensional  (plane  strain). 
A  rigid,  wedge-shaped  indenter  of  angle  0  is  pressed  to  a  depth  8  into  a  calcite  or  sapphire 
substrate  of  nominal  dimensions  50  nm  x  75  nm.  The  indenter  is  rounded  at  the  tip  to  alleviate 
extreme  deformation  that  would  occur  with  a  perfectly  sharp  indenter;  the  analytical  elastic 
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Figure  3.  Wedge  indentation. 


solution  (no  twinning)  for  frictionless  contact  [38^10]  indicates  that  the  pressure  immediately 
under  the  indenter  tip  would  be  infinite  for  a  perfectly  sharp  indenter.  Wedge  angles  of  0  =  90° 
and  0  =  120°  are  considered.  The  shear  strain  under  the  indenter  is  y  ~  tan (90°  —  0 /2),  about 
1  for  the  90°  wedge  and  about  0.58  for  the  120°  wedge.  The  substrate  is  rigidly  bonded  to  the 
indenter  along  the  line  of  contact.  It  is  noted  that  this  boundary  condition  is  an  idealization;  in 
a  real  experiment,  some  slip  would  be  expected  to  occur  between  indenter  and  substrate.  For 
ductile  materials  and  unlubricated  indenters,  a  stick  condition  as  considered  here  may  be  more 
appropriate  than  sliding/frictionless  contact  [40,65].  The  finite  element  software  used  for 
the  simulations  reported  here  does  not  presently  include  the  capability  to  address  multi-body 
contact  or  friction  interactions.  If  and  when  such  capabilities  are  added  in  the  future,  effects 
of  friction  coefficient  on  indentation  can  be  studied. 

Free  boundary  conditions  (i.e.  t  =  0  and  h  —  0)  are  applied  along  the  remainder  of  the 
top  surface  not  in  contact  with  the  wedge.  Recall  from  (20)  that  h  is  the  force  conjugate 
to  variations  in  the  order  parameter;  h  =  0  is  the  logical  boundary  condition  for  a  surface 
free  of  thermodynamic  force  associated  with  the  phase  field,  just  as  t  =  0  is  the  logical 
condition  for  a  surface  free  of  mechanical  force  associated  with  displacements.  These  free 
boundary  conditions  permit  the  order  parameter  and  displacement  to  vary  along  the  surface 
as  an  outcome  of  the  solution.  Condition  h  —  0  is  also  applied  along  the  interface  between 
indenter  and  substrate,  which  permits  formation  of  a  twin  at  the  surface  immediately  beneath 
the  indenter,  as  observed  in  experiments  [1]  (figure  1).  Fixed  boundary  conditions  (i.e.  u  =  0 
and  r]  =  0  on  3  Q)  are  applied  along  the  left,  right  and  bottom  sides  of  the  substrate.  Prescription 
of  alternative  conditions  h  —  0  along  left,  right  and  bottom  sides  did  not  affect  the  solution. 
The  lattice  orientation  of  the  substrate  in  each  case  is  such  that  the  direction  of  twinning  shear 
is  parallel  to  the  direction  of  indentation  loading. 

A  initial  displacement  field  is  applied  to  material  beneath  the  indenter  to  prevent  element 
overlap  or  inversion  during  initial  conjugate  gradient  iterations.  A  defect  ( r\  —  1)  of  small 
size  (radius  ~  5  A)  underneath  the  indenter  is  typically  assigned  as  an  initial  condition.  When 
no  initial  defect  is  prescribed,  twin(s)  do  not  always  nucleate  in  the  numerical  simulations. 
Mathematically,  this  phenomenon  is  understandable  from  examination  of  (22)  and  (24),  which 
can  be  satisfied  when  r\  =  0  everywhere,  even  though  the  total  energy  of  a  system  with  a  twin 
(i.e.  with  r)  >  0  somewhere  in  £2)  may  be  less  than  that  in  which  no  twinning  occurs.  In  other 
words,  condition  r]  =  0  everywhere  may  correspond  to  a  metastable  solution.  The  phenomenon 
is  also  physically  realistic:  Garber  and  Stepina  [8]  found  that  twins  nucleated  immediately  with 
negligible  loads  when  surface  defects  were  present,  but  required  much  greater  indentation  force 
to  nucleate  when  surface  defects  were  absent. 

The  numerical  solution  technique  incorporates  finite  element  discretization  and  conjugate 
gradient  energy  minimization  [20].  Solutions  of  weak  forms  of  equilibrium  equations  in 
section  2  are  obtained,  with  order  parameter  r]  and  displacement  u  the  nodal  degrees  of 
freedom.  For  the  boundary  conditions  prescribed  here,  equilibrium  solutions  corresponds 
to  minima  of  energy  functional  ^  as  in  (42). 
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Table  3.  Phase  field  simulations  (linear). 


Case 

Material 

Twin 

Cabcd 

KAB 

0 

Scale 

1 

CaCC>3 

e+ 

Isotropic 

Isotropic 

90° 

Nanoscale 

2 

CaCC>3 

e+ 

Isotropic 

Isotropic 

120° 

Nanoscale 

3 

CaCC>3 

e+ 

Isotropic 

Anisotropic 

90° 

Nanoscale 

4 

AI2O3 

B 

Isotropic 

Isotropic 

90° 

Nanoscale 

5 

AI2O3 

R 

Isotropic 

Isotropic 

90° 

Nanoscale 

6 

CaCC>3 

e+ 

Isotropic 

Isotropic 

90° 

Lab  scale 

7 

AI2O3 

B 

Isotropic 

Isotropic 

90° 

Lab  scale 

8 

AI2O3 

R 

Isotropic 

Isotropic 

90° 

Lab  scale 

9 

CaCC>3 

e+ 

Anisotropic 

Isotropic 

90° 

Nanoscale 

10 

AI2O3 

B 

Anisotropic 

Isotropic 

90° 

Nanoscale 

11 

AI2O3 

B 

Anisotropic 

Isotropic 

90° 

Lab  scale 

Linear  (i.e.  three  node)  triangular  finite  elements  are  used,  with  significant  mesh  refinement 
underneath  the  indenter.  Element  size  in  such  regions  is  ~  1  A,  sufficient  for  resolution  of 
gradients  in  order  parameter  at  twin-parent  interfaces  by  10+  finite  elements.  Results  of  interest 
are  insensitive  to  further  increases  in  mesh  density.  Indentation  depth  8  is  varied  nominally 
from  0.25  nm  to  4.75  nm  in  increments  of  0.5  nm.  For  each  increment,  initial  conditions  are 
re-applied,  and  then  the  conjugate  gradient  algorithm  is  used  to  determine  equilibrium  order 
parameter  and  elastic  fields. 

Absolute  indentation  force  P  is  computed  by  summing  vertical  nodal  forces  contributing 
to  t2  along  the  contact  line.  Absolute  length  L  of  the  twin  is  determined  by  the  X 2  coordinate 
at  the  tip  demarcated  by  the  local  condition  ri  >  0.1.  In  most  cases,  increasing  the  thickness  of 
the  substrate  above  75  nm  results  in  very  little  change  to  P  or  L;  any  exceptions  are  discussed 
later  in  the  context  of  the  corresponding  results. 

Simulations  discussed  in  sections  4.2  and  4.3  focus  on  nanometer- scale  twins,  i.e. 
nucleation  phenomena.  Length  scales  involved  are  in  agreement  with  previous  studies  of  twin 
nucleation  [2, 20,  36,  37].  Twins  modeled  by  Kosevich  and  colleagues  using  one-dimensional 
theories  [1, 9, 1 1]  are  idealized  with  minimum  thickness  on  the  order  of  the  interplanar  spacing 
(<1  nm),  and  hence  also  describe  nanoscopic  phenomena,  though  the  length  of  twins  studied 
in  such  analytical  models  may  be  orders  of  magnitude  larger.  Using  the  present  phase  field 
approach,  much  larger  (i.e.  laboratory  scale,  dimensions  of  mm)  twins  can  also  modeled  by 
careful  scaling  of  the  size  of  the  indenter  and  substrate,  the  mesh  density  and  the  gradient 
energy  parameter  /c,  such  that  the  finite  element  mesh  is  able  to  resolve  the  thickness  of  twin 
boundary  interfaces.  This  effort  to  model  laboratory  scale  as  opposed  to  nanoscale  twins  is 
described  in  section  4.4. 

Table  3  lists  simulations  discussed  in  forthcoming  parts  of  section  4.  Solutions  obtained 
using  isotropic  elasticity  are  grouped  into  cases  1-8  because  these  will  be  analyzed  together 
later.  Cases  9-11  consider  anisotropic  elasticity. 

4.2.  Isotropic  elasticity 

Contours  of  order  parameter  fields  for  cases  1,  2, 4  and  5  are  shown  in  figure  4.  The  normalized 
indentation  depth  is  A  =  4.75  for  all  cases  in  figure  4.  Normalized  (dimensionless)  depth, 
normalized  indentation  force  per  unit  out-of-plane  length,  and  normalized  twin  length  are 
defined,  respectively,  as 

8  P  L 

A  =  -,  n  =  — ,  A  =  — .  (49) 

l  El  l 
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Figure  4.  Order  parameter  rj  field  for  isotropic  elastic  solutions  at  the  nanoscale  at  indentation 
depth  A  =  4.75:  (a)  case  1:  calcite,  e+  twin,  90°  wedge.  ( b )  Case  2:  calcite,  e+  twin,  120°  wedge, 
(c)  Case  4:  sapphire,  B  twin,  90°  wedge.  ( d )  Case  5:  sapphire,  R  twin,  90°  wedge. 


Effective  Young’s  modulus  E ,  Young’s  modulus  E  and  Poisson’s  ratio  v  are 

-  E  /x(3A.  +  2/r)  X 

E  =  - E  =  - ,  v  =  - . 

1  —  v2  X  +  [i  2(X  +  fi) 


(50) 


The  definition  of  E  is  motivated  by  the  analytical  linear  elastic  solution  (with  no  twinning 
and  frictionless  contact)  [38-40],  wherein  the  mean  contact  pressure  between  indenter  and 
substrate  is  pm  =  E  cot (0/2). 
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Figure  5.  Normalized  indentation  force  versus  indentation  depth  (a)  and  normalized  twin  length 
versus  force  ( b )  for  isotropic  elastic  constants. 


In  each  case,  a  primary  twin  forms  under  the  left  side  of  the  indenter,  where  the  driving 
shear  stress  for  twinning  r  of  (23)  is  positive  in  sign.  Directly  under  the  right  side  of  the  wedge, 
t  is  negative  in  sign,  since  the  shear  strain  is  directed  in  the  anti-twinning  sense  with  respect  to 
the  orientation  of  the  twin  system.  The  sharp  tip  predicted  for  e+  twinning  in  calcite  agrees  with 
observations  from  indentation  experiments  conducted  from  the  1930s  to  the  1970s  [1, 5,  8, 1 1]. 
Long,  lenticular- shaped  micro-twins  in  indentation  of  sapphire  have  been  reported  [32].  The 
predicted  sharp  tip  also  verifies  that  described  by  the  early  analytical  solution  of  Lifshitz  [12]. 

Also  noteworthy  are  the  secondary  twins  that  form  at  the  edge  of  the  indenter  (right  side) 
where  the  shear  deformation  is  directed  in  the  twinning  (as  opposed  to  anti-twinning)  sense. 
This  phenomenon  is  reminiscent  of  the  appearance  of  twins  near  the  outer  edge  of  spherical 
indenters  reported  by  Garber  and  Stepina  [8]. 

Comparing  figure  4(a)  with  figure  4(b),  a  longer  twin  is  produced  in  calcite,  for  the  same 
A,  when  an  indenter  with  larger  wedge  angle  is  used.  Comparing  figure  4(c)  with  figure  4(d), 
for  the  same  A,  a  rhombohedral  (R)  twin  is  significantly  longer  than  a  corresponding  basal 
(B)  twin  in  sapphire.  Recall  from  table  2  that  the  shear  yo  for  B  twinning  is  larger  than  that 
for  R  twinning  (0.635  versus  0.202)  as  is  the  surface  energy  T  (745  versus  125  mJ  m-2).  The 
twinning  shear  for  calcite  is  0.694  (table  1).  Comparing  figure  4(a)  with  figures  4(c)  and  (d), 
for  the  same  A  and  wedge  angle,  e+  twins  in  calcite  and  B  twins  in  sapphire  are  much  closer  in 
length  and  appearance  than  are  R  twins  in  sapphire.  It  is  suggested  that  larger/longer  R  twins 
are  required  to  achieve  comparable  relaxation  in  elastic  strain  energy  than  B  twins  because  of 
the  smaller  stress-free  shear  yo  in  the  former. 

Another  interesting  feature  in  figure  4(b)  is  the  thin  boundary  layer  at  the  surface  of  contact 
under  the  left  side  of  the  wedge  where  r]  <  1 .  This  layer  appears  because  the  imposed  shear 
strain  directly  under  the  120°  indenter  is  about  0.58,  which  is  less  than  the  transformation 
shear  0.694  for  calcite.  On  the  other  hand,  for  indentation  by  90°  wedges,  the  imposed  shear 
strain  is  about  1,  which  is  greater  than  yo,  and  hence  complete  transformation  with  r\  —  1  is 
energetically  favorable  along  the  contact  boundaries  under  the  left  sides  of  the  90°  wedges. 

Indentation  force  versus  depth  is  shown  in  figure  5(a)  for  cases  1-8,  along  with  force 
versus  depth  for  pure  elastic  solutions  (i.e.  no  twinning).  Results  for  no  twinning  are  obtained 
for  90°  wedges.  Twinning  induces  a  significant  reduction  in  force  relative  to  pure  elasticity. 
Use  of  a  120°  wedge  (case  2)  rather  than  a  90°  wedge  produces  a  larger  force  because  of  the 
former’s  larger  contact  area.  In  figure  5(a),  force- versus-depth  data  for  90°  wedge  indentation 
with  twinning  collapse  to  nearly  the  same  curve  when  normalized  to  dimensionless  form  by 
the  appropriate  quantities  via  (49). 
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Twin  length  versus  indentation  force  is  shown  in  figure  5(b).  In  all  cases,  twin 
length  increases  monotonically  with  increasing  force  or  with  increasing  depth  of  indentation. 
Behavior  is  similar  for  e+  twining  in  calcite  and  B  twinning  in  sapphire.  On  the  other  hand,  R 
twins  in  sapphire  are  significantly  longer  under  the  same  normalized  force.  It  was  found  that 
for  twin  lengths  A  <  50,  interaction  with  the  rigid  lower  boundary  dQ  of  the  domain  does  not 
affect  the  twin  length  (recall  r\  =  0  is  imposed  on  X2  =  0).  For  the  R  twin,  the  two  data  points 
corresponding  to  the  largest  values  of  A  indicate  some  influence  of  the  lower  boundary,  which 
tends  to  repel  and  hence  shorten  the  twin.  Perhaps  most  noteworthy  from  figure  5(b),  for  e+ 
twinning  in  calcite  and  B  twinning  in  sapphire,  twin  length  versus  force  for  wedge  indentation 
reduce  to  nearly  the  same  curve  when  normalized  via  (49). 

It  was  remarked  in  section  3.1  that  the  choice  of  anisotropic  surface  energy  in  (43)  might 
facilitate  more  elongated  or  sharper  twins,  as  was  observed  for  the  problem  of  homogeneous 
twin  nucleation  in  previous  work  [20].  Comparing  cases  1  and  3,  use  of  anisotropic  surface 
energy  in  (43)  with  energy  per  unit  area  differing  by  a  factor  of  ~2  in  directions  parallel 
and  perpendicular  to  the  habit  plane  results  in  negligible  differences  in  force  and  length,  and 
contours  of  r\  are  also  negligible  between  cases  1  and  3  (not  shown).  Thus,  the  present 
results  suggest  that  anisotropic  surface  energy  does  not  significantly  affect  predicted  twin 
morphology  during  wedge  indentation,  at  least  for  the  materials  and  ranges  of  physical 
parameters  considered  here. 


4.3.  Anisotropic  elasticity 

Contours  of  order  parameter  rj  are  compared  for  anisotropic  and  isotropic  elastic  constants 
in  figure  6.  Shown  are  results  for  90°  wedges  at  an  indentation  depth  of  A  =  4.75. 
Qualitatively,  with  regards  to  twin  shape,  predictions  of  anisotropic  and  isotropic  models 
agree.  Use  of  anisotropic  elasticity  produces  a  somewhat  shorter  e+  twin  in  calcite  than  does 
the  corresponding  isotropic  model  (figures  6(a)  and  (b)).  Use  of  anisotropic  elasticity  produces 
a  slightly  longer  B  twin  in  sapphire  than  does  the  corresponding  isotropic  model  (figures  6(c) 
and  (d)).  Recall  from  section  3  that  elastic  anisotropy  (in  terms  of  ratios  of  elastic  constants)  is 
significantly  greater  in  calcite  than  sapphire;  however,  relative  influences  of  anisotropy  on  the 
indentation  solutions  might  also  be  influenced  by  other  factors  such  as  differences  in  twinning 
shear  yo  and  different  surface  energies  prescribed  for  the  two  materials. 

Normalized  indentation  force  versus  indentation  depth  for  90°  wedge  indentation  is  shown 
in  figure  7(a)  for  e+  twinning  in  calcite  and  B  twinning  in  sapphire.  All  curves  except  that  for 
anisotropic  calcite  coincide;  the  indentation  force  tends  to  be  lower  in  the  latter.  Normalized 
twin  length  versus  indentation  force  are  shown  in  figure  7(b),  again  for  e+  twinning  in  calcite 
and  B  twinning  in  sapphire.  All  results  collapse  to  nearly  the  same  monotonically  increasing 
(and  nearly  linear)  relationship. 

4.4.  Scaling  to  laboratory  dimensions 

Modeling  twin  nucleation  and  growth  at  dimensions  of  millimeters  or  larger,  i.e.  at  laboratory 
scales  corresponding  to  traditional  indentation  experiments  [1,5,  8, 11],  is  not  feasible  without 
proper  scaling  of  the  interfacial  thickness.  Resolution  of  an  interface  of  thickness  /  of  1  nm, 
as  considered  in  nanoscale  simulations  discussed  in  sections  4.2  and  4.3,  requires  a  numerical 
discretization  with  finite  elements  of  size  on  the  order  of  1  A;  discretization  of  a  domain  of  1  cm2 
would  require  ~1016  such  elements,  for  example.  In  phase  field  calculations,  the  pragmatic 
idea  of  treating  the  interface  as  having  an  adjustable  or  scalable  thickness  has  been  used 
often  for  modeling  larger  systems  [23, 24].  Following  such  an  approach,  here  the  equilibrium 
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Figure  6.  Order  parameter  rj  field  for  isotropic  and  anisotropic  elastic  solutions  at  the  nanoscale 
at  indentation  depth  A  =  4.75  with  90°  wedge:  ( a )  case  1:  calcite,  e+  twin,  isotropic  elasticity. 
(b)  Case  9:  calcite,  e+  twin,  anisotropic  elasticity,  (c)  Case  2:  sapphire,  B  twin,  isotropic  elasticity, 
(i d )  Case  10:  sapphire,  B  twin,  anisotropic  elasticity. 

interfacial  thickness  is  scaled  to  l  =  1  mm,  an  increase  by  a  factor  of  106  above  the  value  of 
1  nm  used  in  sections  4.2  and  4.3.  In  other  words,  in  the  scaled  model,  the  equilibrium  width 
of  the  twin  boundary  is  ~1  mm.  Recall  that  1  nm  corresponds  to  a  distance  over  which  atoms 
deviate  from  their  ideal  positions  in  atomic  simulations  of  twin  boundaries  [22]. 

To  achieve  the  desired  scaling  l  — >  Fl,  values  of  constants  in  (17)  are  adjusted  as 

k  -*  F2k,  a  -*  A,  (51) 
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Figure  7.  Normalized  indentation  force  versus  indentation  depth  (a)  and  normalized  twin  length 
versus  force  ( b )  for  e+  twins  in  calcite  and  B  twins  in  sapphire. 


where  in  this  particular  case,  F  =  106.  The  transformation  in  (51)  preserves  the  ratio  of 
interfacial  energy  to  strain  energy  on  a  per-unit- volume  basis;  W  of  (8)  and  /o  of  (16)  are 
unchanged: 

W  ->  W,  /o  -*  /o.  (52) 

The  gradient  contribution  to  interfacial  energy  per  unit  volume  /  in  (15)  is  also  unchanged 
since 

V^V/F,  K\Vr)\2  ^  (F2K)\(Vr)/n\2  =  K\Vr)\2.  (53) 

In  laboratory  (millimeter)  scale  simulations  discussed  next,  room  temperature  (as  opposed  to 
0  K)  elastic  constants  are  used,  to  better  mimic  ambient  test  conditions. 

Representative  results  are  compared  in  contour  plots  of  r\  in  figure  8.  In  figure  8 (a), 
nanoscale  results  are  shown,  and  absolute  substrate  dimensions  are  50  nm  x  75  nm.  In 
figure  8(b),  laboratory  scale  results  with  parameters  modified  via  (51)  are  shown,  and 
dimensions  are  50  mm  x  75  mm.  Close  similarity  between  cases  1  and  6  in  figures  8(a)  and 
(b)  (calcite,  e+  twin,  isotropic  elasticity)  demonstrates  the  ability  of  the  scaling  to  produce  self¬ 
similar  results  for  linear  elasticity.  Differences  in  twin  length  between  figures  8(a)  and  (b)  are  a 
consequence  of  the  increase  in  compliance  of  the  elastic  constants  at  room  temperature  relative 
to  OK.  As  shown  in  figure  5,  data  for  indentation  force  versus  depth  and  twin  length  versus 
indentation  force  are  nearly  identical  when  normalized  to  dimensionless  form.  For  example, 
compare  results  for  cases  1  and  6  (calcite)  and  cases  4  and  7  (sapphire,  B  twin).  Close  similarity 
between  cases  10  and  11  (sapphire,  B  twin,  anisotropic  elasticity)  demonstrates  the  ability  of 
the  scaling  procedure  to  produce  self-similar  results  for  anisotropic  elasticity. 

Limitations  of  the  above  scaling  method  should  be  noted.  Per  unit  volume,  the  ratio  of 
interfacial  energy  to  strain  energy  is  maintained  by  the  scaling  procedure:  both  W  and  / 
entering  (7)  are  unchanged  by  the  scaling.  This  explains  why  solutions  obtained  by  energy 
minimization  appear  nearly  identical  in  figure  8:  the  same  problem  (with  different  units)  is 
essentially  solved  in  each  case.  However,  per  unit  area,  energy  associated  with  the  interface  is 
increased  by  a  factor  of  F.  Interfacial  energy  per  unit  volume  is  energy  per  unit  area  divided 
by  thickness.  Because  interfacial  thickness  is  increased  by  a  substantial  factor  (i.e.  F ),  small 
features  cannot  be  resolved  in  absolute  dimensions. 

The  present  model  predictions  would  benefit  from  quantitative  comparisons  with 
experiments.  Such  comparisons  are  presently  inhibited  by  limitations  in  reporting  of 
experimental  geometry  and  data.  Regarding  the  former,  it  is  typically  stated  that  a  ‘knife  edge’ , 
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(a)  (b) 


Figure  8.  Order  parameter  rj  field  for  solutions  at  nanometer  and  millimeter  scales  at  indentation 
depth  A  =  4.75  with  90°  wedge:  (a)  case  1:  calcite,  e+  twin,  isotropic  elasticity,  /  =  1  nm. 
( b )  Case  6:  calcite,  e+  twin,  isotropic  elasticity,  /  =  1  mm. 


‘chisel’,  or  ‘wedge’  is  used  to  produce  an  elastic  twin  in  calcite,  but  the  wedge  angle  (i.e.  0 
of  figure  3)  is  not  reported  [1, 9, 7, 11].  Results  of  measurements  of  force  versus  twin  length 
are  given  only  implicitly  as  parameters  entering  dislocation  models  [9, 11];  raw  force  versus 
length  data  enabling  comparison  with  the  present  predictions  are  apparently  unpublished. 
Data  are  available  for  spherical  indentation  of  calcite  [6] ;  however  twinning  under  spherical 
indentation  requires  3D  modeling  outside  the  scope  of  this  work.  New  experiments  studying 
effects  of  wedge  angles  as  reported  in  detail  for  other  elastic  and  plastic  materials  [39]  would  be 
useful.  It  is  reiterated  that  predictions  reported  here  do  agree  qualitatively  with  experimental 
observations:  a  long  thin  twin  forms  under  one  side  of  the  indenter,  the  twin  tip  is  sharp,  and 
the  twin  length  increases  with  increasing  applied  load. 


5.  Phase  field  simulations:  nonlinear  theory 

Simulations  of  indentation  using  the  nonlinear  theory  of  section  2.2  are  described;  numerical 
results  are  analyzed  and  compared  with  those  from  the  linear  theory. 

5.1.  Boundary  value  problem 

The  same  indentation  boundary  value  problems  described  in  section  4  are  considered  here  in  the 
context  of  the  finite  strain  model  of  section  2.2.  The  numerical  solution  technique  incorporates 
procedures  nearly  identical  to  those  discussed  in  section  4.1.  Solutions  of  weak  forms  of  the 
equilibrium  equations  in  section  2.2  are  obtained,  with  order  parameter  r]  and  displacement 
u  =  x  ~  X  the  nodal  degrees  of  freedom.  For  the  boundary  conditions  prescribed  here,  the 
equilibrium  solution  corresponds  to  a  minima  of  energy  functional  as  in  (42).  In  a  few 
exceptional  cases  involving  large  indentation  depths  in  sapphire,  converged  solutions  could 
not  be  obtained  for  all  increments;  in  those  cases,  solution  data  shown  in  subsequent  figures  is 
necessarily  limited  to  that  available  from  converged  solutions.  For  simulations  incorporating 
the  geometrically  nonlinear  neo-Hookean  model,  it  was  found  that  accelerated  convergence 
towards  equilibrium  solutions  could  often  be  obtained  by  using  the  corresponding  linear  elastic 
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Table  4.  Phase  field  simulations  (nonlinear  and  linear). 


Case 

Material 

Twin 

Model 

<P 

Scale 

1 

CaCC>3 

e+ 

Nonlinear 

90° 

Nanoscale 

2 

CaCC>3 

e+ 

Linear 

90° 

Nanoscale 

3 

AI2O3 

B 

Nonlinear 

90° 

Nanoscale 

4 

AI2O3 

B 

Linear 

90° 

Nanoscale 

5 

AI2O3 

R 

Nonlinear 

90° 

Nanoscale 

6 

AI2O3 

R 

Linear 

90° 

Nanoscale 

7 

CaCC>3 

e+ 

Nonlinear 

120° 

Nanoscale 

8 

CaCC>3 

e+ 

Linear 

120° 

Nanoscale 

9 

AI2O3 

B 

Nonlinear 

90° 

Lab  scale 

10 

AI2O3 

B 

Linear 

90° 

Lab  scale 

phase  field  solutions  reported  in  section  4  as  an  initial  guess.  This  initialization  procedure  was 
used  to  obtain  many  of  the  nonlinear  phase  field  solutions  discussed  later  involving  relatively 
large  indentation  depths.  It  was  found  that  the  same  final  equilibrium  solution  was  obtained 
using  this  initialization  method  or  the  initial  conditions  described  in  section  4.1. 

Simulations  discussed  in  sections  5.2  and  5.3  focus  on  nanometer- scale  twins,  i.e. 
nucleation  phenomena,  arising  from  indenters  with  respective  wedge  angles  of  0  =  90°  and 
0  =  120°.  Length  scales  involved  (i.e.  tens  of  nm)  are  similar  to  those  of  prior  analytical  and 
computational  studies  of  twin  nucleation  [2, 20, 36,  37].  Results  of  modeling  laboratory  scale 
as  opposed  to  nanoscale  twins  are  presented  in  section  5.4. 

Table  4  lists  simulations  discussed  in  section  5.  Calcite  (e+  twin  system)  and  sapphire 
(B  and  R  twin  systems)  are  modeled.  Odd-numbered  simulation  cases  incorporate  the 
geometrically  nonlinear  phase  field  model  with  neo-Hookean  elasticity  of  section  2.2.  Even- 
numbered  cases  incorporate  the  geometrically  linear,  linear  elastic  phase  field  model  of 
section  2.1.  For  a  given  nonlinear  solution,  the  linear  solution  for  the  same  material  and 
boundary  conditions  corresponds  to  the  subsequent  case  number  in  table  4. 

5.2.  Nano -indentation  with  90°  wedge 

Contours  of  order  parameter  and  shear  stress  fields  in  calcite  for  simulation  cases  1  (nonlinear 
theory)  and  2  (linear  theory)  are  shown  in  figure  9.  The  normalized  indentation  depth  is 
A  =  4.75.  Normalized  (dimensionless)  depth,  normalized  indentation  force  and  normalized 
twin  length  are  defined,  respectively,  as  in  (49). 

For  each  case  in  figure  9,  a  primary  twin  forms  under  the  left  side  of  the  indenter,  where  the 
driving  shear  stress  for  twinning  r  of  (41)  is  positive  in  sign.  Under  the  right  side  of  the  wedge, 
t  is  negative  in  sign,  since  the  shear  strain  is  directed  in  the  anti-twinning  sense  with  respect 
to  the  orientation  of  the  twin  system.  A  secondary  twin  arises  at  the  edge  of  the  indenter  at  the 
free  surface  (right  side)  where  the  shear  deformation  is  directed  in  the  twinning  (as  opposed 
to  anti-twinning)  sense. 

For  purposes  of  comparison,  contour  legends  for  the  order  parameter  in  these  and 
subsequent  figures  are  restricted  to  r\  e  [0,  1] .  Occasionally,  values  of  r\  slightly  exceeded  unity 
at  one  or  a  few  nodes  immediately  under  the  indenter.  The  current  numerical  implementation  of 
the  phase  field  model  places  no  hard  bounds  on  77;  values  77  <  0  are  set  to  0  in  (5),  while  values 
77  >  1  are  set  to  1  in  (5).  Thus,  values  outside  the  range  [0,1]  offer  no  advantage  with  regards 
to  minimization  of  elastic  strain  energy,  and  are  penalized  in  the  double- well  potential  (16). 

Comparing  figure  9(a)  with  figure  9(b),  order  parameter  contours  and  twin  lengths  are 
similar  for  nonlinear  and  linear  models  in  this  instance.  One  minor  difference  arises:  in  the 
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Figure  9.  Close-up  contours  of  order  parameter  and  shear  stress  for  nano-indentation  with  90° 
wedge  at  depth  A  =  4.75:  ( a )  r],  case  1:  calcite,  e+  twin,  nonlinear  theory.  ( b )  ij,  Case  2:  calcite, 
e+  twin,  linear  theory,  (c)  P12//C  Case  1:  calcite,  e+  twin,  nonlinear  theory.  ( d )  P12//C  Case  2: 
calcite,  e+  twin,  linear  theory. 


nonlinear  solution  (figure  9(a)),  a  thin  layer  of  twinned  material  (r]  ~  1)  emerges  under  the 
right  side  of  the  indenter,  near  the  tip.  Such  a  layer  is  absent  in  the  linear  solution  of  figure  9(b). 
Shear  stress  contours  in  figures  9(c)  and  ( d)  differ  substantially  for  nonlinear  and  linear  models. 
The  shear  stress  is  more  severe,  and  more  strongly  concentrated  under  the  tip  of  the  indenter, 
for  the  nonlinear  model  (figure  9(c))  than  the  linear  model  (figure  9(d)).  Maximum  shear 
stresses  can  be  large,  with  local  magnitudes  exceeding  the  shear  modulus  for  results  of  the 
nonlinear  theory. 
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Figure  10.  Order  parameter  rj  field  for  nano-indentation  of  sapphire  to  depth  of  A  =  3.75  with  90° 
wedge:  (a)  Case  5:  sapphire,  R  twin,  nonlinear  theory  ( b )  Case  6:  sapphire,  R  twin,  linear  theory. 


Contours  of  order  parameter  fields  for  indentation  of  sapphire,  cases  5  (nonlinear,  R  twin) 
and  6  (linear,  R  twin)  are  shown  in  figure  10.  Order  parameter  contours  and  twin  lengths 
are  similar  for  nonlinear  and  linear  models.  In  the  nonlinear  solution,  a  thin  layer  of  twinned 
material  (rj  &  1)  again  emerges  under  the  right  side  of  the  indenter,  near  the  tip,  while  such  a 
layer  is  absent  in  the  corresponding  linear  solution. 

For  the  same  indentation  depth  A,  a  rhombohedral  (R)  twin  is  predicted  to  be  significantly 
longer  than  a  corresponding  basal  (B)  twin  in  sapphire.  Recall  from  table  2  that  the  shear  y0 
for  B  twinning  is  larger  than  that  for  R  twinning  (0.635  versus  0.202)  as  is  the  surface  energy  V 
(745  versus  125  mJ  m-2).  Longer  R  twins  may  be  required  to  achieve  comparable  relaxation  in 
elastic  strain  energy  than  B  twins  because  of  the  smaller  stress-free  shear  yq  in  the  former.  This 
trend  (R  twins  longer  than  B  twins)  was  consistent  among  results  of  both  linear  and  nonlinear 
models. 

Indentation  force  versus  depth  is  shown  in  figure  11(a)  for  cases  1-6.  Force  versus 
depth  data  for  90°  wedge  indentation,  when  normalized  by  the  appropriate  quantities  via  (49), 
are  similar  for  all  cases  at  small  indentation  depths  (A  <  2.5).  At  larger  indentation  depths, 
indentation  forces  for  nonlinear  models  (cases  1 , 3  and  5)  exceed  forces  for  corresponding  linear 
models  (cases  2, 4  and  6).  Such  differences  result,  in  part,  from  the  increase  in  tangent  stiffness 
(e.g.  bulk  modulus)  with  increasing  compressive  pressure  in  the  neo-Hookean  elasticity  model. 
For  e+  twinning  in  calcite  and  B  twinning  in  sapphire,  the  indentation  force  in  the  nonlinear 
model  (case  1  or  3)  only  modestly  exceeds  that  for  the  complementary  linear  model  (case  2  or 
4)  at  A  >  2.5.  The  indentation  force  increases  dramatically  for  case  5  (sapphire,  R  twin)  at 
larger  indentation  depths  relative  to  the  force  for  the  corresponding  linear  model  (case  6).  Order 
parameter  (figure  10)  and  displacement  fields  demonstrate  little  differences  among  linear  and 
nonlinear  solutions  at  large  indentation  depths,  whereas  local  compressive  stress  magnitudes 
become  comparatively  much  larger  for  the  nonlinear  model.  It  is  suggested  that  the  increase  in 
force  in  the  nonlinear  model  is  a  result  of  nonlinear  elastic  contributions  to  stress  (especially 
pressure)  that  become  stronger  at  larger  indentation  depths  and  for  larger  R  twins. 

Twin  length  versus  indentation  force  is  shown  in  figure  11(b).  In  all  cases,  twin 
length  increases  monotonically  with  increasing  force  or  with  increasing  depth  of  indentation. 
Behavior  is  quantitatively  similar  for  e+  twining  in  calcite  and  B  twinning  in  sapphire.  R  twins 
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(a)  Normalized  depth  A  (b)  Normalized  force  n 

Figure  11.  Normalized  indentation  force  versus  indentation  depth  (a)  and  normalized  twin  length 
versus  force  ( b )  for  cases  1-6. 


in  sapphire  are  significantly  longer  under  the  same  normalized  force.  From  figure  11(b),  for 
e+  twinning  in  calcite  and  B  twinning  in  sapphire,  twin  length  versus  indentation  force  data 
normalized  by  (49)  reduce  to  similar  curves. 


5.3.  Nano-indentation  with  120°  wedge 

Contours  of  order  parameter  field  ri  for  cases  7  (calcite,  nonlinear  theory)  and  8  (calcite,  linear 
theory)  for  indentation  with  a  120°  wedge  are  shown  in  figure  12  for  indentation  depths  of 
A  =  3.75  (figures  12 (a)  and  ( b ))  and  A  =  4.75  (figure  12(c)  and  (d)).  At  each  depth,  the 
absolute  length  of  the  long,  thin  primary  twin  is  nearly  identical  for  nonlinear  and  linear  models. 
The  sharp  tip  predicted  for  e+  twinning  in  calcite  agrees  with  observations  from  indentation 
experiments  [1, 5,  8, 11]. 

Predicted  twin  morphologies  are  strikingly  different  for  nonlinear  and  linear  models.  For 
the  linear  theory,  a  single  continuous  primary  twin  forms  under  the  left  side  of  the  indenter.  For 
the  nonlinear  theory,  a  layered  twin  structure  with  alternating  regions  of  large  and  small  values 
of  r\  forms  under  the  left  side  of  the  indenter.  The  number  of  horizontal  layers  increases  with 
increasing  indentation  depth:  e.g.  two  horizontal  layers  are  evident  in  figure  12 (a)  and  three 
arise  in  figure  12(c).  Such  layered  structures  were  not  reported  in  early  experimental  work  on 
indentation  of  calcite  from  the  1930s  to  the  1970s  [1, 5,  8, 1 1];  however,  such  microstructures 
may  simply  have  not  been  detectable  by  the  optical/photographic  equipment  available  at 
that  time.  Kaga  and  Gilman  [10]  observed  layered  (lamellar)  structures  in  twinned  calcite 
specimens  through  etch  pit  studies.  It  was  suggested  that  etch  pit  lines  could  be  associated 
with  positions  where  a  twin  boundary  interface  temporarily  rested  during  growth  of  a  twin  [10]. 
The  alternating  layers  of  twinned  crystal  under  the  indenter  evident  in  the  nonlinear  results 
in  figure  12  are  reminiscent  of  the  finely  twinned  energy-minimizing  structures  in  theories 
of  martensitic  phase  transformations  [27,29].  Because  of  the  energetic  penalty  incurred  by 
boundary  regions  of  r)  /  0,  1  in  the  phase  field  theory,  the  layered  structures  that  emerge  in  the 
nonlinear  theory  must  offer  a  substantial  reduction  in  elastic  strain  energy  relative  to  a  solution 
in  which  a  single  continuous  primary  twin  forms,  as  in  the  linear  theory. 

Other  phase  field  simulations  incorporating  geometrically  linear,  but  possibly  anisotropic, 
elasticity  have  predicted  layered  structures  of  twin  domains  in  the  context  of  martensitic 
transformations  [46, 47].  An  anisotropic  linear  elastic-plastic  finite  element  study  [66]  showed 
that  local  shear  stress  distributions  favor  formation  of  lenticular  shaped  twins,  and  that  above  a 
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Figure  12.  Order  parameter  rj  field  for  nano-indentation  with  120°  wedge:  ( a )  case  7:  calcite,  e+ 
twin,  A  =  3.75,  nonlinear  theory.  ( b )  Case  8:  calcite,  e+  twin,  A  =  3.75,  linear  theory,  (c)  Case  7: 
calcite,  e+  twin,  A  =  4.75,  nonlinear  theory.  ( d )  Case  8:  calcite,  e+  twin,  A  =  4.75,  linear  theory. 


threshold  twinned  volume  fraction,  formation  of  a  layered  microstructure  consisting  of  multiple 
twins  is  energetically  favorable  to  formation  of  a  single  large  twin.  Reported  finite  element 
calculations  did  not  account  for  surface  energy,  which  was  noted  would  become  important  for 
very  small  twins  in  the  nucleation  stage  [66] .  None  of  these  other  works  addressed  indentation 
loading. 

In  all  cases  shown  in  figure  12,  a  thin  boundary  layer  arises  at  the  surface  of  contact  under 
the  left  side  of  the  wedge,  where  r]  <  1 .  This  layer  appears  because  the  imposed  shear  strain 
directly  under  the  120°  indenter  is  about  0.58,  which  is  less  than  the  transformation  shear 
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Figure  13.  Close-up  contours  of  order  parameter  and  shear  stress  for  nano-indentation  with  120° 
wedge  at  depth  A  =  4.75:  ( a )  r],  case  7:  calcite,  e+  twin,  nonlinear  theory.  ( b )  r],  Case  8:  calcite, 
e+  twin,  linear  theory,  (c)  Pn/li,  Case  7:  calcite,  e+  twin,  nonlinear  theory.  ( d )  P12//X,  Case  8: 
calcite,  e+  twin,  linear  theory. 


0.694  for  calcite.  On  the  other  hand,  for  indentation  by  90°  wedges  discussed  in  section  5.2 
(figure  9),  the  imposed  shear  strain  is  about  1,  which  is  greater  than  yo,  and  hence  complete 
transformation  with  r]  =  1  is  energetically  favorable  along  the  contact  boundaries  under  the 
left  sides  of  90°  wedges.  Because  the  layered  twin  structure  of  the  nonlinear  theory  does 
not  emerge  in  simulations  of  indentation  by  90°  wedges,  the  appearance  of  such  layers  for 
indentation  by  120°  wedges  (imposed  shear  of  0.58  <  yo  =  0.694)  could  be  explained  by 
similar  arguments. 

Close-up  views  of  order  parameter  r]  and  shear  stress  component  Pn  are  shown  in  figure  13 
for  cases  7  and  8.  Shear  stress  contours  in  figures  13(c)  and  (J)  differ  substantially  for  nonlinear 
and  linear  models.  As  was  observed  in  figure  9,  the  shear  stress  is  more  severe,  and  more 
strongly  concentrated  under  the  tip  of  the  indenter,  for  the  nonlinear  model  (figure  13(c))  than 
the  linear  model  (figure  13 (J)).  Alternating  layers  of  shear  stress  that  accompany  the  lamellar 
twin  structure  predicted  by  the  nonlinear  theory  are  just  visible  in  figure  13(c). 

Indentation  force  versus  depth  is  shown  in  figure  14 (a)  for  cases  7-10.  Twin  length  versus 
indentation  force  is  shown  in  figure  14 (b).  Results  for  cases  7  and  8  are  discussed  here;  cases 
9  and  10  are  discussed  in  section  5.4.  Force-versus-depth  data  for  120°  wedge  indentation, 
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o  1  2  3  4  5  0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4 

(a)  Normalized  depth  A  (b)  Normalized  force  n 

Figure  14.  Normalized  indentation  force  versus  indentation  depth  (a)  and  normalized  twin  length 
versus  force  ( b )  for  cases  3,  4  and  7-10. 


when  normalized  by  (49),  are  almost  identical  for  linear  and  nonlinear  models.  Indentation 
force  and  twin  length  predicted  by  nonlinear  theory  exceed  corresponding  values  predicted  by 
linear  theory  by  less  than  1%.  Twin  length  increases  monotonically  with  increasing  force  or 
with  increasing  depth  of  indentation. 


5.4.  Scaling  to  laboratory  dimensions 

For  scaling  according  to  the  transformation  /  ->  Fl,  values  of  constants  in  (17)  are  adjusted 
as  in  (51)  of  section  4.4.  Again,  for  the  present  problem,  F  =  106,  scaling  the  interface 
thickness  from  1  nm  to  1  mm.  In  laboratory  (millimeter)  scale  simulations  discussed  next, 
room  temperature  (as  opposed  to  0  K)  elastic  constants  are  used,  to  better  represent  ambient 
experimental  conditions  [1, 5,  8, 11]. 

Representative  results  are  compared  for  nonlinear  and  linear  theories  in  contour  plots  of 
T)  in  figure  15,  corresponding  to  basal  twins  in  sapphire  induced  by  90°  wedge  indentation. 
Figures  15 (a)  and  ( b )  can  be  compared  with  nanoscale  results  in  figures  6(c)  and  ( d ),  though 
the  indentation  depth  is  slightly  smaller  in  the  former.  In  figures  6(c)  and  (d),  nanoscale  results 
are  shown,  and  the  absolute  substrate  dimensions  are  50  nm  x  75  nm.  In  figures  15(a)  and 
(b),  laboratory  scale  results  with  parameters  scaled  via  (51)  are  shown,  and  the  dimensions 
are  50  mm  x  75  mm.  Note  that  twin  morphologies  are  similar  for  model  predictions  at  each 
length  scale.  Linear  and  nonlinear  theories  also  predict  similar  order  parameter  profiles  in 
figure  15.  Similar  to  what  was  observed  in  figure  9(a),  only  in  the  nonlinear  solution  does 
a  thin  layer  of  twinned  material  (rj  ~  1)  emerge  under  the  right  side  of  the  indenter  near 
the  tip. 

Indentation  force  versus  depth  is  shown  in  figure  14 (a);  twin  length  versus  indentation 
force  is  shown  in  figure  14(b).  Force- versus-depth  data  for  90°  wedge  indentation  of  sapphire 
with  B  twins,  when  normalized  by  the  appropriate  quantities  via  (49),  is  quantitatively  very 
similar  for  linear  and  nonlinear  models  at  both  the  nanoscale  and  laboratory  scale  (cases  3,  4, 
9,  and  10).  Normalized  twin  length  data  for  linear  and  nonlinear  models  at  each  length  scale 
also  collapse  to  nearly  the  same  monotonically  increasing  curve,  as  demonstrated  by  the  nearly 
indistinguishable  results  for  sapphire  in  figure  14(b).  Close  similarity  among  results  for  cases 
3  and  9  in  figure  14  demonstrates  the  ability  of  the  scaling  method  of  section  4.4  to  produce 
geometrically  self-similar  results  for  nm-  and  mm-sized  specimens  using  the  nonlinear  theory. 
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Figure  15.  Order  parameter  rj  field  for  laboratory  (mm)  scale  indentation  with  90°  wedge  at  depth 
A  =  4.25:  (a)  case  9:  sapphire,  B  twin,  nonlinear  theory.  ( b )  Case  10:  sapphire,  B  twin,  linear 
theory. 


6.  Conclusions 

A  phase  field  model  for  twinning  in  elastic  crystals  has  been  developed  and  exercised  for  the 
problem  of  indentation  with  wedges.  Two  transparent  materials  have  been  studied:  calcite  and 
sapphire.  Predictions  of  nonlinear  (i.e.  finite  strain)  theory  have  been  compared  with  those  of 
geometrically  linear  theory. 

Phase  field  predictions  of  indentation  agree  qualitatively  with  experimental  observations:  a 
long  thin  twin  forms  asymmetrically  beneath  one  side  of  the  indenter,  the  tip  of  the  twin  is  sharp 
and  the  length  of  the  twin  increases  with  increasing  force.  Normalized,  dimensionless  variables 
for  indentation  depth,  indentation  force,  and  twin  length  have  been  derived.  It  is  emphasized 
that  physically  realistic  twin  shapes  are  predicted  by  a  model  whose  only  parameters  are  the 
elastic  constants,  twin  boundary  energy,  and  twin  boundary  thickness,  all  of  which  can  be 
obtained  from  independent  experiments  or  quantum/molecular  mechanics  calculations. 

Qualitatively  similar  twin  shapes  have  been  obtained  using  isotropic  and  anisotropic  elastic 
constants.  The  difference  in  predicted  twin  length  between  isotropic  and  anisotropic  models  is 
greater  in  sapphire  than  in  calcite.  Basal  and  rhombohedral  twins  have  been  studied  in  sapphire: 
predicted  rhombohedral  twins  are  longer  than  basal  twins  for  the  same  indentation  force.  Use 
of  anisotropic  rather  than  isotropic  surface  energy  has  little  effect  on  twin  morphology  or 
indentation  force  for  the  particular  indentation  boundary  conditions  and  ranges  of  material 
parameters  considered  here. 

A  scaling  method  has  been  developed  for  modeling  behavior  of  specimens  of  (mm-to- 
cm)  sizes  on  the  order  of  those  studied  experimentally  in  traditional  indentation,  as  opposed 
to  nano-indentation.  Nearly  identical  twin  morphologies  are  obtained  for  nanometer-scale 
specimens  and  millimeter  (lab)  scale  specimens.  Twin  length  versus  indentation  force  data, 
when  properly  normalized  to  dimensionless  form,  collapse  to  nearly  the  same  monotonically 
increasing  curve  for  e+  twins  in  calcite  and  basal  twins  in  sapphire,  for  nm-  and  mm-scale 
simulations,  and  for  90°  and  120°  wedges. 

Indentation  forces  are  greater  in  the  nonlinear  model  than  the  linear  model  because  of 
the  decreasing  elastic  compliance  with  increasing  pressure  in  nonlinear  model.  Normalized 
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relationships  between  twin  length  and  indentation  force  are  similar  for  linear  and  nonlinear 
theories  at  both  nanometer  and  millimeter  scales.  Twin  morphologies  are  similar  for  linear  and 
nonlinear  theories  for  indentation  with  90°  wedges:  in  each  case,  a  single,  continuous  primary 
twin  forms  under  one  side  of  the  indenter,  and  a  small  secondary  twin  forms  at  the  free  surface 
at  the  opposite  edge.  Perhaps  most  interestingly,  in  the  nonlinear  model,  indentation  of  calcite 
with  a  120°  wedge  produces  a  lamellar  twin  structure  between  the  indenter  and  the  long  sharp 
primary  twin.  The  number  of  twin  lamellae  increases  with  increasing  indentation  depth.  This 
complex,  layered  microstructure  is  not  predicted  by  the  linear  theory,  which  instead  predicts  a 
single  continuous  primary  twin. 
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